In [ ]:
import numpy as np
import pandas as pd
from tqdm import tqdm
import plotly.graph_objects as go
import plotly.express as px
import plotly.io as pio
pio.renderers.default = "notebook+pdf"
System Variables
In [ ]:
# Event codes
ARRIVAL=0
DEPARTURE=1
event_codes = {ARRIVAL: "ARRIVAL", DEPARTURE: "DEPARTURE"}
# Patient status codes
PENDING=0
ADMITTED=1
RELOCATED=2
REJECTED=3
DISCHARGED=4
status_codes = {PENDING: "PENDING", ADMITTED: "ADMITTED", RELOCATED: "RELOCATED", REJECTED: "REJECTED"}
# Simulation parameters
SIM_TIME = 31 # Simulation time in days
# Ward parameters
WARDS = dict(
names = ['A', 'B', 'C', 'D', 'E', 'F'],
lams = [14.5, 11.0, 8.0, 6.5, 5.0, 13.0],
mu_invs = [2.9, 4.0, 4.5, 1.4, 3.9, 2.2],
urgency_points = [7, 5, 2, 10, 5, 0],
)
NUM_WARDS = len(WARDS['names'])
RELOCATION_PROBABILITIES = np.array([
[0.0, 0.05, 0.10, 0.05, 0.80, 0.00],
[0.2, 0, 0.50, 0.15, 0.15, 0.00],
[0.30, 0.20, 0, 0.20, 0.30, 0.00],
[0.35, 0.30, 0.05, 0, 0.3, 0.00],
[0.20, 0.10, 0.60 ,0.10, 0, 0.00],
[0.20, 0.20, 0.20, 0.20, 0.20 ,0]
])
Penalty Computation
In [ ]:
def penalty(urgency_points, counts):
return np.dot((counts[RELOCATED] + counts[REJECTED]), urgency_points)
Simulation Setup
In [ ]:
def simulation_setup(sim_time=SIM_TIME):
# Patients
patients = dict( ID=[], type=[], ward=[], status=[], U1=[], U2=[], burn_in=[] )
# Precompute arrival and departure times for each ward
events = dict( event_type=[], patient_ID=[], time=[] )
patient_id = 0
for i, (lam, mu_inv) in enumerate(zip(WARDS['lams'], WARDS['mu_invs'])):
# Sample patients for ward
clock = 0
while clock < sim_time:
# Sample arrival time
U1, U2 = np.random.uniform(size=2)
clock += -np.log(U1)/lam
# Add patient to event list if arrival time is before simulation end
if clock <= sim_time:
# Event
events['time'] += [clock, clock-np.log(U2)*mu_inv]
events['patient_ID'] += [patient_id]*2
events['event_type'] += [ARRIVAL, DEPARTURE]
# Patient
patients['ID'].append(patient_id)
patients['type'].append(i)
patients['ward'].append(i)
patients['status'].append(PENDING)
patients['U1'].append(U1)
patients['U2'].append(U2)
patients['burn_in'].append(0)
# Update patient ID
patient_id += 1
# Sort events by time
idx = np.argsort(events['time'])
for key in events.keys():
events[key] = [events[key][i] for i in idx]
return patients, events
Inspect the Event List
In [ ]:
# Setup simulation
patients, events = simulation_setup()
# Create dataframes and map codes to names
patient_df = pd.DataFrame(patients); patient_df['status'] = patient_df['status'].map(status_codes)
patient_df['type'] = patient_df['type'].map({i: WARDS['names'][i] for i in range(NUM_WARDS)})
event_df = pd.DataFrame(events); event_df['event_type'] = event_df['event_type'].map(event_codes)
# Print dataframes
print("Patients:")
print(patient_df.head().to_string(index=False), "\n")
print("Events:")
print(event_df.head().to_string(), "\n")
Patients: ID type ward status U1 U2 burn_in 0 A 0 PENDING 0.675212 0.036522 0 1 A 0 PENDING 0.914004 0.612931 0 2 A 0 PENDING 0.268128 0.037591 0 3 A 0 PENDING 0.582789 0.018506 0 4 A 0 PENDING 0.414830 0.484123 0 Events: event_type patient_ID time 0 ARRIVAL 1011 0.017203 1 ARRIVAL 0 0.027085 2 ARRIVAL 451 0.029163 3 ARRIVAL 1 0.033286 4 ARRIVAL 765 0.038836
Simulation Loop
In [ ]:
def simulate(patients, events, capacities, burn_in_period=15):
# Clear patient status
N = len(patients['ID'])
patients['ward'] = [None]*N
patients['status'] = [PENDING]*N
# List for saving ward states
states = []
# Dict for saving counts
counts = {
ADMITTED: np.zeros(NUM_WARDS, dtype=int),
RELOCATED: np.zeros(NUM_WARDS, dtype=int),
REJECTED: np.zeros(NUM_WARDS, dtype=int),
}
# Simulate loop
occupancy = np.zeros(NUM_WARDS, dtype=int)
for i, (event_type, patient_id) in enumerate(zip(events['event_type'], events['patient_ID'])):
# Check if patient is in burn-in period
if events['time'][i] < burn_in_period:
patients['burn_in'][patient_id] = True
# Get patient type
patient_type = patients['type'][patient_id]
if event_type == ARRIVAL:
# Admit patient
if occupancy[patient_type] < capacities[patient_type]:
ward = patients['type'][patient_id]
status = ADMITTED
else:
# Relocate patient (or reject if alternative ward is full)
ward = np.random.choice(NUM_WARDS, p=RELOCATION_PROBABILITIES[patient_type])
status = RELOCATED if occupancy[ward] < capacities[ward] else REJECTED
# Update patient
if status != REJECTED:
occupancy[ward] += 1
patients['ward'][patient_id] = ward
patients['status'][patient_id] = status
if events['time'][i] >= burn_in_period:
counts[status][patient_type] += 1
# Discharge patient
elif patients['status'][patient_id] != REJECTED:
occupancy[patients['ward'][patient_id]] -= 1
states.append(occupancy.copy())
events['states'] = states
# Check if simulation ended with non-zero occupancy
if any(occupancy != 0):
print('Error: Simulation ended with non-zero occupancy')
return dict(states=states, counts=counts, penalty=penalty(WARDS['urgency_points'], counts))
Run Simulation
In [ ]:
patients, events = simulation_setup()
capacities = [49, 28, 22, 16, 16, 34]
sim_out = simulate(patients, events, capacities)
states = sim_out['states']
fig = go.Figure()
for i, state in enumerate(np.array(events['states']).T):
fig.add_trace(go.Scatter(x=events['time'], y=state, mode='lines', name=WARDS['names'][i]))
fig.update_layout(title='Ward occupancies', xaxis_title='Time', yaxis_title='Occupancy', width=800)
fig.show()
Gradient Computation
In [ ]:
def compute_gradients(capacities):
M = len(capacities)
grads = np.zeros(M)
n = 10
def simulate_penalty(capacity_adjustment):
total_penalty = 0
for _ in range(n):
capacities[i] += capacity_adjustment
patients, events = simulation_setup()
sim_out = simulate(patients, events, capacities)
total_penalty += sim_out['penalty']
capacities[i] -= capacity_adjustment
return total_penalty / n
for i in range(M):
p_small = simulate_penalty(-1)
p_large = simulate_penalty(1)
grads[i] = (p_large - p_small) / 2
return grads
def gradient_descent(capacities, beds_to_F=34):
penalties = np.zeros(beds_to_F+1)
F_rel_prob = np.zeros(beds_to_F+1)
# Use the samee simulations to estimate performance metrics
sim_setup_samples = [simulation_setup() for _ in range(100)]
# Realocate beds to F
for i in tqdm(range(beds_to_F+1)):
# Simulate and compute performance metrics
sim_out = [simulate(patients, events, capacities) for patients, events in sim_setup_samples]
penalties[i] = np.mean(np.array([out['penalty'] for out in sim_out]))
num = [out['counts'][RELOCATED][-1] + out['counts'][REJECTED][-1] for out in sim_out]
den = [out['counts'][ADMITTED][-1] + out['counts'][RELOCATED][-1] + out['counts'][REJECTED][-1] for out in sim_out]
F_rel_prob[i] = np.mean([num/den for num, den in zip(num, den)])
# Compute gradients and update capacities
grads = compute_gradients(capacities)
idx = np.argmax(grads[:-1])
capacities[idx] -= 1
capacities[-1] += 1
# Correct for last iteration
capacities[idx] += 1; capacities[-1] -= 1
return dict(capacities=capacities, penalties=penalties, F_rel_prob=F_rel_prob)
Optimizing Bed Allocation
In [ ]:
# Run gradient descent for different bed totals (Takes a few minutes)
# out165 = gradient_descent([55, 40, 30, 20, 20, 0])
# capacities165 = out165['capacities']
# out170 = gradient_descent([55, 40, 30, 20, 20, 5])
# capacities170 = out170['capacities']
# out180 = gradient_descent([55, 40, 30, 20, 20, 15])
# capacities180 = out180['capacities']
# Found capacities
capacities165 = [49, 28, 22, 16, 16, 34]
capacities170 = [48, 31, 25, 18, 14, 34]
capacities180 = [52, 35, 25, 16, 18, 34]
In [ ]:
out165 = {'capacities': [50, 32, 25, 14, 10, 34],
'penalties': np.array([ 994.7 , 992.51, 1023.78, 1016.06, 1031.28, 1054.56, 1037.07,
1060.79, 1059.14, 1096.48, 1074.49, 1084.01, 1090.24, 1089.09,
1106.09, 1110.83, 1113.56, 1110.67, 1121.34, 1120.74, 1126.26,
1139.59, 1153.75, 1174.76, 1190.52, 1155.77, 1203.88, 1227.05,
1269.71, 1258.28, 1284.74, 1305.39, 1327.39, 1327.1 , 1349.44]),
'F_rel_prob': np.array([1. , 0.96497114, 0.93127537, 0.89639315, 0.8651794 ,
0.83104137, 0.7953772 , 0.76118428, 0.72714368, 0.69326257,
0.66094461, 0.62798531, 0.59402437, 0.5632216 , 0.53226308,
0.49888117, 0.46506489, 0.43620332, 0.40787952, 0.37578122,
0.34587529, 0.31573004, 0.2874799 , 0.2606288 , 0.23393766,
0.20971284, 0.18552666, 0.16272165, 0.14025691, 0.11959395,
0.09943001, 0.0831262 , 0.06695472, 0.05308433, 0.04153909])}
Theoretical Optimal Number of Beds in F
In [ ]:
def erlangB_formula(m):
if m == 0:
return 1
lam = 13; s = 2.2
A = lam*s
def power_div_factorial(i):
return np.exp(i*np.log(A) - np.sum(np.log(range(1, i+1))))
return power_div_factorial(m)/(np.sum([power_div_factorial(i) for i in range(m+1)]))
In [ ]:
fig = go.Figure()
fig.add_trace(go.Scatter(x=np.arange(35), y=[erlangB_formula(m) for m in range(35)], mode='lines', name='Erlang B'))
fig.add_trace(go.Scatter(x=np.arange(35), y=out165['F_rel_prob'], mode='lines', name='Simulation'))
fig.update_layout(title='Probability of relocating in ward F', xaxis_title='Number of beds in ward F', yaxis_title='Probability', width=800)
fig.show()
fig = go.Figure()
fig.add_trace(go.Scatter(x=np.arange(35), y=out165['penalties'], mode='lines', name='Simulation'))
fig.update_layout(title='Penalty as a function of number of beds in ward F', xaxis_title='Number of beds in ward F', yaxis_title='Penalty', width=800)
fig.show()
Estimate Rates for Admission, Relocation and Rejection
In [ ]:
def estimate_rates(samples, capacities):
num_samples = len(samples)
rates = np.zeros((num_samples, 3, 6))
for i in range(num_samples):
patients, events = samples[i]
sim_out = simulate(patients, events, capacities)
counts = sim_out['counts']
rates[i, 0] = counts[ADMITTED]
rates[i, 1] = counts[RELOCATED]
rates[i, 2] = counts[REJECTED]
rates /= rates.sum(axis=(1, 2))[:, None, None]
rates = rates.mean(axis=0)
rates = pd.DataFrame(rates.T, index=WARDS['names'], columns=['ADMITTED', 'RELOCATED', 'REJECTED'])
return rates
In [ ]:
num_samples = 100
sim_setup_samples = [simulation_setup() for _ in range(num_samples)]
for n, cap in zip([165, 170, 180], [capacities165, capacities170, capacities180]):
rates = estimate_rates(sim_setup_samples, cap)
print(f'{n} beds:')
print(rates.to_string(), '\n')
165 beds: ADMITTED RELOCATED REJECTED A 0.206927 0.016383 0.029591 B 0.101809 0.043963 0.041178 C 0.056730 0.051762 0.029614 D 0.085913 0.014185 0.010681 E 0.039231 0.022740 0.022262 F 0.216036 0.006077 0.004917 170 beds: ADMITTED RELOCATED REJECTED A 0.208288 0.015399 0.029214 B 0.114944 0.040273 0.031734 C 0.065890 0.045475 0.026742 D 0.097142 0.008047 0.005590 E 0.036114 0.026882 0.021237 F 0.216036 0.006548 0.004446 180 beds: ADMITTED RELOCATED REJECTED A 0.225689 0.011821 0.015391 B 0.127633 0.034712 0.024605 C 0.070692 0.048506 0.018909 D 0.094211 0.010523 0.006044 E 0.049021 0.020295 0.014916 F 0.216036 0.007356 0.003638